The Projected Gross-Pitaevskii Equation for harmonically confined Bose gases 
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We extend the Projected Gross Pitaevskii equation formalism of Davis et al. [Phys. Rev. Lett. 87, 160402 
(2001)] to the experimentally relevant case of harmonic potentials. We outline a robust and accurate numerical 
scheme that can efficiently simulate this system. We apply this method to investigate the equilibrium properties 
of a harmonically trapped three-dimensional Bose gas at finite temperature, and consider the dependence of 
condensate fraction, position and momentum distributions, and density fluctuations on temperature. We apply 
the scheme to simulate an evaporative cooling process in which the preferential removal of high energy particles 
leads to the growth of a Bose-Einstein condensate. We show that a condensate fraction can be inferred during 
the dynamics even in this non-equilibrium situation. 

PACS numbers: 03.75.-b, 03.75.Hh 



I. INTRODUCTION 

The near zero temperature formalism for Bose-Einstein 
condensates is based upon the well-established Gross- 
Pitaevskii equation (GPE) that describes the macroscopically 
occupied condensate orbital. It has been proposed that the 
GPE can also be used to model the non-equilibrium dynam- 
ics of finite temperature Bose gases Ql yl 0. 01 ■ The essen- 
tial idea is that highly occupied modes of a quantum field 
are well described by a classical field, as is well-known for 
the case of electromagnetic fields. Recently, calculations 
have been performed by several groups using classical fields 
B S II Ull 03 [H Qj and have shown the usefulness of 
this approach. One of the key advantages of the classical field 
method is that it is non-perturbative; however care must be 
taken to appropriately cutoff the spectrum to avoid an ultravi- 
olet catastrophe, analogous to that occurring in the Rayleigh- 
Jeans theory of blackbody radiation. Ideally this cutoff can be 
made from a priori thermodynamic analysis of the system. In 
the finite temperature formalism of Davis et al. I6l 171 fl3ll . the 
cutoff is explicitly incorporated through the use of a projec- 
tor that is diagonal in the single particle basis of the system 
Hamiltonian. This ensures that a consistent energy cutoff is 
established, and provides a natural separation of the system 
into regions of low energy modes (the highly occupied co- 
herent region) and high energy modes (the sparsely occupied 
incoherent region). Ignoring the incoherent region and its cou- 
pling to the coherent region, the equation of motion for the low 
energy modes is termed the Projected Gross-Pitaevskii Equa- 
tion (PGPE). By itself the PGPE can only provide a partial 
description of the system, however it contains the modes that 
are most significantly modified by the effects of interactions 
and are difficult to include quantitatively in traditional kinetic 
theories. Thus the PGPE by itself can be a useful tool for 
providing insight into the evolution of ultra-cold Bose gases. 

The Stochastic GPE formalism developed by Gardiner and 
coworkers 1 14, 15] proposes a practical scheme for coupling 
the coherent and incoherent regions in a manner suitable for 
non-equilibrium calculations. To account for the interaction 
with the incoherent region, this formalism introduces noise 



terms into the PGPE that transfer energy and particles into 
and out of the coherent region. 

A related technique originating in quantum optics is the 
phase space method known as the truncated Wigner approx- 
imation. This was first applied to ultra-cold Bose gases by 
Steel et al. [1 1611 . and further developed by Sinatra and co- 
workers 

E3 uE 0- 

A number of recent calculations using 
this method can be found in Refs. |20, 2lL l22ll23ll24ll . 

The main purpose of this paper is to develop the PGPE ap- 
proach for the experimentally relevant case of three dimen- 
sional harmonic traps. To do this we have adapted a recent 
numerical scheme by Dion and Cances Il25ll by introducing 
an explicit energy cutoff in the single particle (harmonic os- 
cillator) basis. We use the scheme to examine the proper- 
ties of trapped gases as a function of the temperature, such 
as condensate fraction and density fluctuations. While some 
of these properties were studied in Ref. [10], we believe that 
our numerical method is superior and our results are free from 
the significant uncertainties found in this paper. Finally, as a 
non-equilibrium application we use the PGPE to simulate the 
formation of a Bose-Einstein condensate in an evaporatively 
cooled thermal cloud. 



II. FORMALISM 

The theoretical formalism we numerically implement in 
this paper has been developed in Refs. 0, [3. IT3l . and 
here we briefly summarise the main points of this formalism, 
adapted for application to inhomogeneous systems. 

A dilute Bose gas is well-described by the second quantized 
Hamiltonian 



H — H sp + Hi, 



(1) 



where 



Hi = 



dx$t(i) ( -i—V 2 + V trap (x) ) V(x), (2) 
2m 



2 



are the single particle and interaction Hamiltonians respec- 
tively, and ip(~x.) is the quantum Bose field operator that anni- 
hilates a particle at position x. The inclusion of an external 
trapping potential, assumed to be harmonic in form i.e. 

Wx) = \mu? z {X 2 x x 2 + X 2 y f + S 2 ), (4) 

is the major new feature of this paper as compared to the ear- 
lier numerical investigations made by one of us in Refs. |00- 
In Eq. (|3} particle interactions have been approximated by a 
contact interaction of strength Uq = 4Trh 2 a/m, where m is 
the atomic mass, and a is the s-wave scattering length. The 
harmonic trap geometry is defined by the angular oscillation 
frequencies along each axis, i.e. lo x , u y , and u> z . For conve- 
nience we express the x and y frequencies relative to u> z , by 
introducing the relative frequency parameters X x = uo x /uo z 
and X y = LUy/uj z . 

The exact equation of motion for the field operator is given 
by the Heisenberg equation of motion 



in — - 

dt 



+ VWp(x) ) V(x)+«7 ^ t (x)^(x)^(x). 



(5) 

For situations of experimental relevance the Hilbert space is 
enormously large, and directly solving Eq. (0 is not possible 
without further approximation. The essence of our approach 
is to split the field operator into two parts representing the 
coherent and incoherent regions. We define the projection op- 
erators 



P{F(i)} = 5> n (x) f d 3 xV:(x')f(x'), (6) 
Q{F(Si)} = £>„(x) f rf 3 xV;(x')^(xO, (7) 

where n G C defines the modes that make up the coherent 
region C, and ip n (x) is the nth eigenfunction of the basis that 
diagonalizes single particle Hamiltonian H sp . We define 



*(x) = F{^(x)}, 
ry(x) - QMx)}, 



(8) 



which we refer to as the coherent (^) and incoherent (?)) field 
operators respectively. This division is based on the average 
occupation of the states. The coherent field "J is chosen to 
describe the low-lying highly occupied states of the system, 
i.e. modes containing of order five or more particles. The in- 
coherent field 7/(x) contains the complementary states which 
are sparsely occupied. Our particular interest is in situations 
near equilibrium, where this separation can be conveniently 
introduced via thermodynamic arguments by making an ap- 
propriate energy cutoff E cut in the single particle spectrum. 

To derive the PGPE, we apply the projection operator (|6j to 
the equation of motion for the field operator l|5}- The funda- 
mental approximation we make, often referred to as the classi- 
cal field approximation, is to neglect the quantum mechanical 



nature of coherent field operator, i.e. set ^(x) — ► Vf^x) (a c- 
number function) due to the high occupation numbers of these 
modes. By making this approximation, the equation of motion 
(|5} is transformed to a form involving couplings between the 
coherent and incoherent fields of some complexity, (e.g. see 
lll3lo . As a further approximation we neglect the interaction 
between the coherent and incoherent regions and simply con- 
sider the equation of motion for (x) in isolation 



ih 



a*(x) 

dt 



-^V 2 + Krap(x) j *!Xi 



Pi t/ |*(x)| 2 *(x) 



(9) 



which is the Projected Gross-Pitaevskii equation (PGPE). 
This equation is of a similar form as the usual Gross-Pitaevskii 
equation, however the classical field 'I'(x) has a distinct in- 
terpretation: It represents the quantum field for many low- 
lying modes, rather than just the condensate mode. Several 
substantial approximations have been made to reduce the full 
Heisenberg equation of motion for the field l|5} to the PGPE 
(|9}, however previous studies have demonstrated that this final 
form contains a rich set of physics (e.g. see 



in. NUMERICAL APPROACH 

The modes of the system are of central important in the as- 
sumptions used to derive the PGPE, and care must be taken in 
numerical implementations to ensure the modes are faithfully 
represented. It is in our opinion that any useful simulation 
technique must satisfy the following requirements: 

1. The space spanned by the modes of the simulation 
should match that of the coherent region of the phys- 
ical system being simulated as closely as possible. That 
is, the modes should be the single particle modes of the 
system up to the prescribed energy cutoff E cut . 

2. The assumption of high occupancy in all modes neces- 
sitates that the numerical scheme must propagate all 
modes accurately. 

Most commonly used methods for propagating Schrodinger- 
type equations do not satisfy these requirements, in particular 
many methods do not propagate all modes of the numerical 
basis faithfully. This leads to negligible errors if the highest 
modes are unoccupied, as is the case for the T = GPE. 
However, it is clear that methods based on such assumptions 
will not be appropriate for simulating the PGPE. 

Before introducing the numerical scheme used in thispa- 
per, we review how the method used by Davis et al. (a IH 
to simulate homogeneous Bose gases addresses the aforemen- 
tioned conditions. For the homogeneous system the modes are 
plane-wave like and are suitable to grid (Fourier) methods of 
propagation (also used in 1 12]). To define the coherent region 
Davis et al. instigated an energy cutoff by using an explicit 
projection operator in momentum space, ensuring condition 1 
was satisfied. To satisfy condition 2, sufficiently many states 
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outside the cutoff should be retained to ensure that the non- 
linear terms in the evolution equation were exactly evaluated 
(without aliasing) for all modes of the coherent region. This 
was not realized at the time and aliasing is present in the calcu- 
lations presented in |fj \Ji |8|] . However, the presence of alias- 
ing mainly modifies the action of the nonlinear term and does 
not affect the equilibrium state or the conclusions reached in 
these papers. For further details we refer the reader to |26]. 

We comment that simply using an unmodified grid method 
as in IT2I1 does introduce a cutoff into the system in momen- 
tum space, but does not satisfy the two criterion listed above. 
First, the energy cutoff is anisotropic, varying in magnitude by 
a factor of three with direction in momentum space (in three 
dimensions). Second, the largest momentum states will be 
aliased in the calculation of nonlinear terms and their dynam- 
ics will be misrepresented. The situation is even worse when 
using an unmodified grid method to simulate classical field 
dynamics for a trapped Bose gas as in fiol Hill . It is this 
issue that we address here. 



A. Brief review of numerical method 



with the value of _E cut chosen to be appropriate for the phys- 
ical system under consideration. For later convenience we 
define n cu t to be the largest index occurring in C, i.e. the 
quantum number of the highest energy oscillator state in the 
coherent region (i.e. n cut « E cut /huj z ). 

In the basis representation, the PGPE JlQi takes the form 



dc h 



dt 



+ CFi 

nut (*)] , (14) 



where 



F lmn {^) = J d 3 x^(x)^(yK(z)|vI/(x ; i)| 2 vI/(x,t), 

(15) 

is the matrix element of the nonlinear term. An important 
observation made in Ref. 12511 is that these matrix elements 
(11 5i can be computed exactly with an appropriately chosen 
Gauss-Hermite quadrature. To show this we note that be- 
cause the harmonic oscillator states are of the form (p n (x) ~ 
H n (x) exp(— x 2 /4), where H n (x) is a Hermite polynomial 
of degree n, the wave function can be written as 



The method we have used to simulate the Projected Gross- 
Pitaevskii equation with an explicit cutoff in en ergy de- 
rives from a recent numerical scheme by Dion et al. 12511 . We 
briefly review our adaption of this method. 

To begin we rescale the Projected Gross-Pi taevskii equation 
(|9} by introducing units of distance xq = - l[ /H/2mu> z , time 
to = a;" 1 , and hence energy hw z , with uj z the trap frequency 
along the ^-direction. With these choices, we have 



9* 



1 



(A 2 z 2 + Xly 2 + z 2 )* + C|#| 2 f JO) 



where we have defined the nonlinear coefficient as C = 
NcUo/hui z XQ, and for clarity have used untilded variables to 
indicate quantities expressed in computational units. We take 
the wave function to be normalized to unity, so that the total 
number of atoms within the coherent region, Nc, appears in 
the definition of the nonlinearity constant. To simplify our dis- 
cussion of the numerical method, we will take the harmonic 
trapping potential to be isotropic, i.e. A^ = X y = 1. This 
allows us to avoid using cumbersome notation to account for 
different spectral bases in each direction. 
The classical field ^'(x, t) is expanded as 

(t) <pi(x)(p m (y)<Pn(z), (li) 

{l,m,n}£C 

where {(p n (x)} are the eigenstates of the ID harmonic oscil- 
lator Hamiltonian satisfying 



dx 2 + A X 



ip n {x) = e n <p n (x), 



(12) 



with eigenvalue e n = (n + 3). The energy cutoff is imple- 
mented by restricting the summation indices to the set 



C = {l,m,n : Ei + e m + e n < E cut }, 



(13) 



*(x) = Q{x,y,z)e 



-(x 2 +y 2 + Z 2 )/4 



(16) 



where Q(x, y, z) is a polynomial that, as a result of the cutoff, 
is of maximum degree n cu t in the independent variables. 

It follows that because the interaction term d!5> is fourth 
order in the wave function, it can be written in the form 



<f xe -(* 2 +y 2 +* 2 ) 



P{x,y,z) 



(17) 



where P(x, y, z) is a polynomial of maximum degree 4 n cut 
in the independent variables. Identifying the exponential 
term as the usual weight function for Gauss-Hermite quadra- 
ture, the integral can be exactly evaluated using a three- 
dimensional spatial quadrature grid of 8 (n C ut) 3 points. Thus 
we have verified that the matrix elements i 7 i mn (2} can be ex- 
actly calculated. We refer the reader to Ref. 12511 for more 
details of how to efficiently implement the spatial transforma- 
tion and numerical quadrature. 



IV. SIMULATION PROCEDURE 

A. Micro-canonical ergodic evolution 

The evolution of the classical field preserves several con- 
stants of motion. These can be considered as macroscopic 
parameters that constrain the micro-states available to the sys- 
tem. For the PGPE, the most important such constant of mo- 
tion is the total energy, given by the energy functional 



- / d 3 x **(x) ( -^V 2 + y trap (x) ) *(x) 



-*7 |*(x 



-\|4 



(18) 
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Another important constant of motion is the field normaliza- 
tion, given by 

iV[#] = / d 3 i \y(x)\ 2 . 



As discussed in section HIl Al we take the classical field to be 
normalized to unity for the initial condition: a scaling choice 
that causes the coefficient of the nonlinear term in the PGPE 
to be proportional to the initial number of particles in the co- 
herent region. The PGPE may have other constants of motion, 
such as angular momentum components in traps with the ap- 
propriate symmetry, however here we will only consider sit- 
uations where these are approximately zero and can be ne- 
glected. 

The lowest energy solution to the energy functional is given 
by the time-independent Gross-Pitaevskii equation (e.g. see 



Ref. I27ln . We denote the energy of this solution E g . This 
solution corresponds to T = 0, however this situation lies 
outside the validity regime of the PGPE since only a single 
mode, i.e. the condensate mode, is highly occupied. For ap- 
plication of the PGPE our interest is in regimes with E > E g 
such that the system is at finite temperature with many highly 
occupied modes. 




Ex" 3 ] 
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Figure 1 : Thermalized classical fields for two values of energy. Den- 
sity slices taken in the y = plane for a system in an anisotropic trap 
with X x — V& and X y = 1. (a) Low energy case with an average 
energy per particle of E = 10hui z . (b) Higher energy case with an 
average energy per particle of E — 24hui z . For both simulations 
C = 2000 and the ground state energy is E g = 8.54/k<j z . The single 
particle cutoff energy is E cut = 31hui z , below which 1739 single 
particle states remain in the coherent region. 

For the simulations we present here we make use of the 
ergodic hypothesis, which we discuss further below. How- 
ever, an immediate consequence of ergodicity for the study 



of equilibrium properties is that the precise details of the ini- 
tial conditions in a simulation are irrelevant. In practise we 
choose initial conditions to provide the desired values for the 
constants of motion. The amplitudes for each single parti- 
cle mode have a random phase and occupation, but are con- 
strained to fix the overall normalization to unity and the en- 
ergy to the desired value. In general such an initial choice 
will not be a typical equilibrium state, but under evolution the 
system rapidly thermalizes. This thermalization process has 
been investigated for the homogeneous case in Q. In Fig. ^ 
we show typical density profiles of thermalized classical fields 
for the harmonically trapped system. The cases in Fig. [Qa) 
and (b) differ in the energy of the fields. It is clear from these 
figures that the classical field has a somewhat chaotic appear- 
ance, and this undergoes constant evolution as the constituent 
modes mix through the nonlinear interaction. 



B. Time-averaging correlation functions 

The energy functional dl 81 is nonlinear, and it is not fea- 
sible to determine the entire phase space of classical field 
configurations consistent with a particular choice of energy. 
This prohibits the calculation of quantities using an ensemble 
averaging approach. However, we can determine correlation 
functions for the system by using the ergodic hypothesis, i.e. 
replacing ensemble averages by time-averages following the 
prescription 



ensemble 



lim 



drF c [*(x,t)] 



where ^[^(x)] is some functional (correlation function) of 
the field. 

The general spatial correlation functions we are interested 
in are usually expressed as an ensemble average over a product 
of quantum field operators such as 



F Q = (*T( Xl ) . . . *T(x*)*(x J+1 ) . . . *(x„)>. 

Note that we have only considered correlation functions that 
involve the field operator for the coherent region. Also we 
have restricted our attention to same-time correlations, how- 
ever in principle multi-time correlation functions could also 
be computed. To evaluate these correlation functions within 
the framework of the classical field theory we make the substi- 
tution ^(x) — ► ^(x), as discussed in section |n] and replace 
ensemble averaging with time averaging. The resulting ex- 
pression that we evaluate numerically is 



(F c [*(x)]) timeBV0 . = -L^FcM^tj 

where {tj} is a set of N s equally spaced time instances at 
which the classical field has been calculated. For this choice to 
well-approximate the ensemble average we require N s 3> 1, 
and the time span over which averaging is done to be long 
compared to the slowest time scale in the problem, e.g. the 
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longest harmonic oscillator period. For the equilibrium results 
we present in this paper these conditions are well satisfied: we 
use 1400 discrete samples of the classical field taken over an 
evolution time of approximately 191 oscillator periods (i.e. 
t Ns -t x = \2W/u z ). 

C. Temperature 

For comparison with experiments and other theories it is 
crucial to be able to identify the temperature of the classical 
field simulations rather characterize a particular result by its 
energy. Previous attempts to determine temperature have been 
based on fitting the occupation of high energy modes to pertur- 
bative calculations for the spectrum based on Hartree-Fock- 
Bogoliubov (HFB) theory (3 0] . For harmonically trapped 
gases, calculation of the HFB modes is much more difficult, 
and limits temperature calculations to perturbative regimes. 
Goral et al. 11(11 have estimated the temperature in harmon- 
ically trapped classical field simulations in a manner analo- 
gous to that done in experiments, by fitting the high momen- 
tum components of the system to a noninteracting distribution. 
The results of that analysis suffered from excessively large er- 
rors, and indicate this approach would not be useful for any 
systematic investigation of thermodynamic properties. In ad- 
dition, the high energy modes of these calculations are un- 
likely to have been represented accurately 12811 . 

In recent work 1 8] one of us has generalized Rugh's dynam- 
ical definition of temperature |29] to the PGPE. This scheme 
has the advantage that it is non-perturbative, and is quite accu- 
rate. This scheme can be extended to the harmonically trapped 
case and is used to calculate the temperature of the simula- 
tions presented in this paper. Because the implementation of 
this scheme in the harmonically trapped case is a trivial exten- 
sion to the homogeneous implementation, we refer the reader 
to Ref. 1 8] for details. 

V. RESULTS 

For the results presented in this paper we have simulated 
a 3D system with = v°, X y = 1 and an energy cutoff 
of _E cut = 3 lfhu z . This choice leads to a coherent region 
containing 1739 harmonic oscillator modes. 

A. Condensation 

1 . Condensate fraction 

Identification of the condensate fraction and mode function 
for an inhomogeneous system with interactions is nontrivial. 
This issue was addressed by Penrose and Onsager in 1956 
I3CB . who extended the concept of Bose-Einstein condensa- 
tion from an ideal gas to the case of superfluid helium. Their 
primary criterion for condensation is that a single eigenvalue, 
no, of the one-body density matrix, p(x, x'), becomes an ex- 
tensive parameter of the system. With regard to obtaining a 



quantitative description of the condensate, they also showed 
that no and its corresponding eigenvector are the condensate 
occupation and mode respectively. 

The one-body density matrix can be written as an ensemble 
average of the quantum field operators as follows: 

pfox') = (^(X^x'^cnscmblc (19) 

As we neglect the incoherent region in this paper, we restrict 
our consideration to the one-body density matrix for the co- 
herent region, i.e. /9 C (x,x') = (* t (x)£(x')) cnscm bic, and 
using the procedure outlined in section lTV Bl we calculate this 
matrix as 

1 N s 

pc(x,x')« — 5>*(x,t,)*(x',^). (20) 

S 3=1 

It is more convenient numerically to represent the density ma- 
trix in the spectral representation, 



Pij = c i (* n ) c i (*" ) ' (2 1 ) 

s n—1 

where Cj(t n ) are the spectral amplitudes of the classical field 
at time t n , and the index j labels all three quantum numbers 
needed to specify the oscillator mode that Cj refers to. The 
efficiency of the spectral representation affords us the ability 
to work with the entire one-body density matrix. The one- 
body density matrix was also time-averaged in the grid based 
method reported in Ref. 11(11 . however their analysis was lim- 
ited to the s-wave component. 




0.002 0.004 0.006 0.008 0.01 0.012 
Temperature [T Q ] 

Figure 2: Equilibrium condensate occupation as a function of tem- 
perature for a range of interaction strengths. The C — case is 
computed using the equipartition relation for the modes of the co- 
herent region only. The dashed lines are simple polynomial fits to 
the numerical results. Other simulation parameters are the same as 
in Fig. □ The averaging was performed using N s = 1400 samples 
over a evolution period of T — 1200/cJz . The unit of temperature is 
To — Nhui z /kB Oili where fcs is the Boltzmann constant. 
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In Fig. |2 we show the condensate fraction (f c = n^/Nc) 
as a function of system temperature for a range of interaction 
strengths. The rather distinctive shape of the curves arises 
because these calculations are for fixed Nc and E cut . For a 
real system with a fixed total number of particles, the portion 
of atoms in the coherent region and the cutoff energy used 
to define the coherent region will change with temperature. 
These issues will be important considerations in making ex- 
perimental comparisons, however in this paper we concern 
ourselves with characterizing properties of the classical field 
method and will only consider the case of Nq and E cut as 
being fixed and independent of temperature. 

For the interacting homogeneous Bose gas there is no mean 
field shift of the critical temperature. Instead critical fluctua- 
tions give the leading order shift with interaction strength to 
the transition temperature, and this shift is upwards. In the 
trapped system there is a mean field effect on the transition 
which reduces the critical temperature 13111 . In Fig. [5]we ob- 
serve this downward shift in the transition temperature as the 
interaction strength is increased, in qualitative agreement with 
recent experiments by Gerbier et al. [32]. We refer the reader 
to [33] for the application of the PGPE as a quantitative model 
of these experiments. 

a. Suitable averaging to determine the condensate frac- 
tion: Using linear algebra arguments it can be shown that 
the condensate fraction determined according to the prescrip- 
tion we have outlined above will have a lower bound of 

f c >max{l/N„l/Gch ( 22 ) 

where N s is the number of samples used to construct the den- 
sity matrix in Eq. (1211 . and Gc is the number of single particle 
states in the coherent region. Typically we take N s < Gc so 
that 1 /N s forms a lower bound for the condensate fraction. 
Equality in the bound holds if the individual classical fields 
included in the time-average are mutually orthogonal, i.e if 
J2i c *i (tm)cj(t n ) = 6 mn V{m,n : 1 < m,n < N s }. Re- 
sult (1221 implies, for instance, that to determine a condensate 
fraction below 1% will require us to take N s > 100 samples. 



2. Condensation: influence on density distributions in momentum 
and position space 

It is interesting to consider how the presence of a con- 
densate affects the position and momentum density pro- 
files of the system. Indeed, it was the appearance of an 
anisotropic peak in the momentum distribution that was used 
as one of the signatures of condensation in the first exper- 
iments 1 34]. In Figs. |3ja)-(c) we compute the momen- 
tum column density for cases above and below the transi- 
tion temperature. In detail, the quantity we calculate is the 
momentum space column density given by J dk y |$(k, t)\ 2 
(i.e. integrated along the fcj,-direction), where $(k, t) — 
(27r)~ 3 / 2 J d 3 x exp(— ik ■ x)4 r (x, t) is the momentum space 
wave function and k = (k x , k y ,k z ). Time-averaging this, i.e 



calculating 




yields the average column density shown in Figs. |3Ja)-(c). 
The peak momentum density of the three cases considered 
varies over a wide range and for presentation clarity we have 
used a logarithmic density scale in Figs. |3Ja)-(c). The appear- 
ance of a narrow peak in the momentum distribution for the 
condensed state is clearly observed in Figs. |3jb) and (c), how- 
ever the logarithmic density scale somewhat suppresses the 
prominence of this feature: the peak momentum column den- 
sity in Fig. |5Ic) is 26 times larger than the peak density in Fig. 
12a). We also note that the momentum distribution changes 
from being isotropic in Fig. [3Ja) where there is no condensate 
(as calculated according to the criterion in Sec. IV A 11 1 to ex- 
hibiting distinctive anisotropy for the condensate momentum 
peak in Figs. |3jb) and (c). This anisotropy is directly related 
to the ratio of the trap frequencies. The usual experimental 
method for imaging the momentum distribution is to take an 
absorption image of the system after allowing it freely expand 
in the absence of the trap. Interaction effects during expansion 
significantly suppresses the contrast in momentum widths of 
condensed and uncondensed systems. However, the contrast 
has been revealed by experiments using Bragg spectroscopy 
techniques [35] that are able to probe the momentum distribu- 
tion in situ. 

Similarly we can construct the column density distribution 
in position space as 




This is equivalently obtained by integrating out the y-direction 
of the diagonal one-body density matrix pc{ x , x ) (1201 . These 
distributions are shown (beneath the associated momentum 
distribution) in Figs. 0d)-(f). These results emphasize 
that while the momentum distribution undergoes substantial 
changes at the transition as discussed above, the position dis- 
tribution changes in a much more subtle manner. The calcu- 
lations presented in 1 10] looked only at the behavior of slices 
through the position distribution of the field. 



B. Density Fluctuations 

Determining the condensate fraction using the Penrose- 
Onsager criterion is a probe of first-order coherence in the 
system. Indeed, the existence of a condensate is equivalent 
to off-diagonal long-range order, i.e. the system is spatially 
coherent. To fully characterize the field it is necessary to con- 
sider higher-order correlations in the system. 

To demonstrate the usefulness of the PGPE, we use it to 
calculate the normalized n-th order coherence function at zero 
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Figure 3: Time-averaged column densities in momentum (a)-(c) and position (d)-(e) space of classical field simulations with various energies. 
Cases: (a) and (d) E = 24hoj z ; (b) and (e) E = 18fioj z ; (c) and (f) E = 10huj z . Other simulation parameters are the same as in Fig. □ The 
averaging was performed using N s — 1400 samples over a evolution period of T — 1200/oj z . 



spatial separation, defined as 



3n(x) 



£ f (x)) (*(x 



(*t( x )*(x))' 



(23) 



Once again we note that we have restricted our attention to 
correlation functions involving the coherent field operator. 
The normalized coherence functions have been calculated for 
the case of n = 2 by Dodd et al. |36] using Hartree-Fock- 
Bogoliubov theory in the Popov approximation, which should 
be valid for large condensate fractions. In contrast the classi- 
cal field result should be applicable as long as our assumption 
of high mode occupancy is satisfied. 

Using the ergodic averaging procedure (see section I1VB> 
we evaluate d23l for the cases n = 2 and n = 3. The results, 
shown in Fig. |4j are for the case of C = 2000 in the pancake 
geometry trapping potential (\ x = v° and X y = 1). For ref- 
erence, we note that Fig. Q(a) corresponds to a single profile 
used for the results in Figs. |4](a) and (d), and similarly Fig. ^ 
(b) corresponds to a single profile used for the results in Figs. 



0Jc) and (f). These coherence functions were evaluated by an- 
gular averaging about the symmetry axis in the x — plane 
(in addition to the time-averaging), and are plotted against the 
radial distance from the trap center in that plane. 

The coherence functions provide a useful characterisation 
of quantum fields. For instance g n = 1 for coherent fields 
such as a laser, whereas g n = n\ for thermal light fields. These 
features are clearly apparent in our results for the matter-wave 
field. When a condensate is present it occupies the center of 
the trapping potential and dominates the thermal fraction of 
atoms in this region. This is clearly seen in Figs. |4](a) and 
(d), where for a condensate fraction of 87%, there is a shape 
transition from coherent behavior (g n = 1) near the trap cen- 
ter, to thermal behavior (<?„ = n!) at a radius of approximately 
r = 5xq. 

In Figs. |4jb) and (e) the same qualitative behavior is seen, 
however the smaller condensate fraction (approximately 24% 
for this case) causes the boundary between coherent and ther- 
mal regions to be closer to the trap center. Also, the coher- 
ence near the trap center is suppressed (i.e. slightly increased 
from unity), indicating that the dominant thermal fraction of 
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Figure 4: Normalized coherence functions Q2 (r) and gs, (r) (see text) 
as a function of radial position in the x — plane (averaged over 
angle about the symmetry axis), and for various condensate fractions. 
The simulation parameters are the same as in Fig. Q and all results 
are for the case C = 2000. For reference the coherent value (g„ — 1) 
and thermal value (g n = n\) of the coherence functions are indicated 
as dotted and dashed lines in the plots respectively. 
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the system is penetrating into the region of the condensate, or 
is in some manner causing increased fluctuations in the con- 
densate mode. 

The results shown in Figs. 0](c) and (f) are for sufficiently 
high energy that the condensate fraction is zero. For this case 
we find the thermal behavior (g n — nl) is present everywhere 
in the system. 



Figure 5: Evolution of an evaporatively cooled matter- wave field. 
(a)-(d) Momentum space density in the k y = plane, (e)-(h) corre- 
sponding position space density in the y = plane. The + signs are 
used to indicate thezero-coordinate in the plots for reference. The 
initial state for the simulation is shown in Fig. 0b), and the evapora- 
tion is applied until t — 64tt/ui z by setting *&(x) = for |z| > 9xo 
during the simulation. Other parameters are as in Fig.Q 



A. Evaporative cooling simulation 



VI. APPLICATION TO NON-EQUILIBRIUM DYNAMICS 

The most interesting application of the PGPE technique 
will be to the non-equilibrium Bose gas. In many cases this 
will require a quantitative description of the incoherent region, 
which is beyond the scope of this paper, but is the subject of 
on going work 1 15]. To give a qualitative demonstration of a 
non-equilibrium application, we consider a simplified simu- 
lation of evaporative cooling using the PGPE. We emphasise 
that this is not intended for quantitative comparison with ex- 
periment, but the calculation highlights several interesting fea- 
tures of the non-equilibrium dynamics, particularly in relation 
to the identification of the condensate. 



To perform a simulation of evaporative cooling, we begin 
with the classical field in an equilibrium state above the tran- 
sition temperature. The cooling is implemented in a man- 
ner analogous to that used in experiments: high-energy atoms 
that are able to venture into regions far from the trap center 
are selectively removed. We do this by absorbing the por- 
tion of the classical field which extends outside the spatial 
region \z\ < 9xq (i.e. setting it to zero at each time step 
of the simulation). We note that for the initial state consid- 
ered [which is the same as the state shown in Fig. QJb)] only 
a small fraction of the field extends into this region, and so 
the normalisation and energy of the classical field is lost rel- 
atively slowly during the cooling [also see Fig. |6jb)]. After 
32 trap periods the cooling mechanism is turned off and the 
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system is allowed to evolve and rethermalize for a subsequent 
20 trap periods. A summary of the classical field dynamics at 
instances during this simulation is shown in Fig. [5] The initial 
momentum and position profiles are shown in Figs. [2a) and 
(e) respectively. After approximately 20 trap periods of cool- 
ing a momentum peak has developed in the distribution near 
k = [see Fig. |3b)]. During these early stages of growth the 
condensate undergoes strong sloshing and breathing dynam- 
ics as fierce mixing occurs between the forming condensate 
and other low-lying quasiparticle modes. The images in Figs. 
|5jc) and (g) show the field at the end of the evaporative cool- 
ing (32 trap periods). These figures show a large condensate 
centered about zero momentum [see Fig. EJc)] and a rela- 
tively settled position distribution [see Fig. |5jg)]. The con- 
densate exhibits breathing dynamics, however this is signifi- 
cantly quenched relative to the strong dynamics seen at earlier 
stages of condensate growth. The condensate at the end of the 
rethermalization period is shown in Figs. |5Jd) and (h). 



B. Time-dependent condensate fraction 

>From the momentum space images of Fig. |U it seems 
quite obvious when a Bose-Einstein condensate has formed. 
A sharp peak suddenly appears in momentum space, whereas 
there is no such clear signature in the real space distributions. 
We note that the momentum space images are on a logarith- 
mic scale — so the peak is even more obvious using a lin- 
ear scale. However, in Ref. [10] it was stated that the time- 
averaging inherent in the imaging process of real experiments 
was essential for the splitting of the system into a condensed 
and non-condensed fraction. Our results are in clear disagree- 
ment with this conclusion, and it is at least qualitatively ap- 
parent that condensation has occurred from a single image of 
the classical field. 

To quantitatively investigate this observation and to exam- 
ine the growth of the condensate, we first apply the Penrose- 
Onsager approach discussed in section IV A II to the evapora- 
tive cooling simulation. The cooling is only carried out in one 
dimension, and so the dynamics proceed rather slowly. Thus it 
seems that we should be able to estimate the one-body density 
matrix at a given time by time-averaging over short periods. 
We calculate the condensate fraction at trap period intervals by 
averaging the one-body density matrix over that interval (by 
summing the classical field at 30 discrete instances during that 
interval), and diagonalizing it. The results are shown in Fig. 
|6ja) as the open circles. Due to the finite time over which we 
are able to average, the initial condensate fraction calculated is 
non-zero despite the initial state having zero condensate frac- 
tion. Because the system is not in equilibrium during the evap- 
oration it is not clear that the Penrose-Onsager approach is ap- 
plicable, however, the characteristic s-shaped curve we find in 
Fig. |6ja) is expected theoretically and has been observed ex- 
perimentally (see ll37ll and references therein). Since the evap- 
orative cooling mechanism is dissipative, particles and energy 
are lost from the system. The evolution of these quantities in 
the simulation are shown in Fig. |6jb). This shows that during 
the cooling 74% of the particles and 88% of the energy in the 
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Figure 6: Growth of condensate and loss of energy and normaliza- 
tion during evaporative cooling process, (a) The condensate frac- 
tion of the total remaining classical field. The open circles are cal- 
culated by diagonalizing the one-body density matrix estimated by 
time-averaging over two trap periods, while the solid curve is cal- 
culated by fitting bimodal distributions to single-shots as described 
in the text, (b) The classical field energy (dashed line) and normal- 
ization (solid line) as a function of time. The point of time at which 
evaporation is stopped is marked by a vertical dotted line. Simulation 
parameters are explained in Fig. [5] 



classical field is lost. 

The second approach we take is to fit bimodal distributions 
to single-shot column densities of the classical field in mo- 
mentum space. This computational method mimics the ac- 
tual experimental procedure that is used for fitting condensate 
plus thermal cloud absorption images in the laboratory. The 
fitting function is the sum of two Gaussian profiles of differ- 
ing widths, and three separate least-squares fits are carried out 
along each of the x, y, and z axes. Our fitting procedure de- 
termines the boundary of the condensate, and the condensate 
number is the integral of the column density in this region less 
the estimated density of the thermal cloud. It seems beneficial 
to fit column densities, as this averages overs some of the fluc- 
tuations apparent in slices through planes of the classical field, 
in a manner reminiscent of the spatial averaging carried out in 
the realisation of a single trajectory in 12311 . A visualisation of 
the results of this fitting method are displayed in Fig. 

The results of this second method are also shown in Fig. 
|6ja) as the solid curve, and the results are in remarkable agree- 
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Figure 7: Example of the bimodal fitting procedure for a single-shot 
of a classical field with a condensate fraction of about 0.24. (a) The 
image on the left is the column density of the classical field, while 
the image on the right is the bimodal fit. (b) A slice through another 
column density of the same field. The open circles are the data points 
and the solid line is the fitted curve. 



ment with the Penrose-Onsager approach. Thus it seems to us 
that the condensate fraction in a static harmonic trap can be es- 
timated from a single-shot image of the classical field, without 
any time-averaging being necessary. In hindsight this seems 
obvious, as this is the standard experimental procedure and it 
seems to have had some success. The time-averaging that oc- 
curs due to the finite exposure time in experimental imaging 
is only on the order of tens of microseconds 13811 - which is 
almost instantaneous on the time-scale of the matter-wave dy- 
namics. In particular this must be the case for non-destructive 
techniques such as phase-contrast imaging otherwise little in- 
formation would be gained. Intuitively the Penrose-Onsager 
approach should require averaging over a time-scale of order 
a trap period. Thus it seems to us that the claim that time- 
averaging is necessary to identify the condensate in Ref. lUoll 



is incorrect. Indeed, if imaging was performed using longer 
duration exposures, photon recoil effects on the atoms would 
dominate over any averaging of the matter wave dynamics 
(e.g. see Sec. 3.5.2 of Ref. 113911 1. A question of interest is 
how the methods we have used here would perform in strongly 
non-equilibrium situations, however we leave this investiga- 
tion for future work. 

As a final remark, we again wish to emphasise that the sim- 
ulation example neglects physical processes which would be 
important in a quantitative model of an evaporative cooling 
experiment. This arises through our ignorance of the inco- 
herent region which would be responsible for the significant 
transfer of particles and energy into the coherent region. How- 
ever, this model does help to illustrate the rather complex 
dynamics that occur in the coherent region as it responds to 
the selective removal of high energy components. The PGPE 
method models the complete non-perturbative dynamics of 
the low-lying modes and is naturally suited to considering 
non-equilibrium situations such as evaporative cooling. 



VII. CONCLUSIONS 

In this paper we have presented an efficient numeri- 
cal scheme for implementing the Projected Gross-Pitaevskii 
equation formalism in three-dimensional harmonic traps with- 
out any axis of symmetry. The main feature of this scheme is 
that it implements a consistent energy cutoff in the harmonic 
oscillator basis and is suitable to efficient and accurate nu- 
merical simulation on modern computer workstations. As an 
application of the method we have used it to simulate a fi- 
nite temperature Bose gas in an anisotropic harmonic trap both 
above and below the critical temperature. Using the ergodic 
hypothesis we have obtained equilibrium quantities such as 
the condensate fraction and the temperature for these simula- 
tions, and have calculated the second- and third-order normal- 
ized coherence functions. As a non-equilibrium application 
we have used the PGPE to simulate the growth of a conden- 
sate from an evaporatively cooled thermal cloud. We have 
managed to identify the condensate fraction in this calcula- 
tion from both the diagonalization of the time-averaged den- 
sity matrix, as well as single-shot column densities in momen- 
tum space of the classical field. 
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